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Abstract 

Complex non-linear relationships exist between air and soil temperature responses to climate change. Despite its influence 
on hydrological and biogeochemical processes, soil temperature has received less attention in climate impact studies. Here 
we present and apply an empirical soil temperature model to four forest sites along a climatic gradient of Sweden. Future 
air and soil temperature were projected using an ensemble of regional climate models. Annual average air and soil 
temperatures were projected to increase, but complex dynamics were projected on a seasonal scale. Future changes in 
winter soil temperature were strongly dependent on projected snow cover. At the northernmost site, winter soil 
temperatures changed very little due to insulating effects of snow cover but southern sites with little or no snow cover 
showed the largest projected winter soil warming. Projected soil warming was greatest in the spring (up to 4°C) in the 
north, suggesting earlier snowmelt, extension of growing season length and possible northward shifts in the boreal biome. 
This showed that the projected effects of climate change on soil temperature in snow dominated regions are complex and 
general assumptions of future soil temperature responses to climate change based on air temperature alone are inadequate 
and should be avoided in boreal regions. 



Citation: Jungqvist G, Oni SK, Teutschbein C, Putter MN (2014) Effect of Climate Change on Soil Temperature in Swedish Boreal Forests. PLoS ONE 9(4): e93957. 
doi:10.1371/journal.pone.0093957 

Editor: Shibu Jose, University of Missouri, United States of America 

Received October 19, 2013; Accepted March 11, 2014; Published April 18, 2014 

Copyright: © 2014 Jungqvist et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits 
unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. 

Funding: This study was funded by Swedish science council, Formas, SKB, Kempe foundation, Mistra and SLU. This study is part of larger forest water programs; 
ForWater and Future Forest, studying the effect of climate and forest managements on forest water quality. The ENSEMBLES data used in this work was funded by 
the EU FP6 Integrated Project ENSEMBLES {Contract 505539) whose support is gratefully acknowledged. MNF was supported by the ECCO and DomQua projects. 
The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. 

Competing interests: The authors have declared that no competing interests exist. 

* E-mail: stephen.oni@slu.se 



Introduction 

There is increasing consensus that the global climate is getting 
warmer in comparison to pre-industrialized times [1,2]. This is 
supported by a growing number of observations; changing weather 
patterns [3], increasing water temperatures [4,5], rising sea levels 
[6], melting of ice and snow in the arctic or subarctic regions [7,8] 
amidst other signals. The effects of global warming are not 
restricted to air temperature alone but can also change precipi- 
tation patterns as well as soil temperature. Snow is a very effective 
insulator and snow cover can effectively decouple air and soil 
temperature during the winter. Thus, soil temperature response to 
climate change may differ from that of air temperature depending 
on changes in the timing and duration of the winter snowpack [9] . 

Despite their importance in controlling watershed biogeochem- 
ical and hydrological processes, soil temperature data are generally 
less available than air temperature measurements. Soil tempera- 
ture controls biogeochemical processes such as dissolved organic 
carbon export [10], length of growing season [11,12], rates of 
mineralization [13,14] or decomposition of soil organic matter 
[15] and nutrient assimilation by plants [16,17], weathering of 
base cations [18] as well as forest productivity [19]. 

Soil temperature is influenced by many factors. The insulating 
effects of snow as well as changes in the timing and intensity of 
snowfall can have significant feedback effects on soil temperature 
dynamics in the northern boreal landscape during winter and 



spring [20] . Studies have shown the possibility of colder soils in a 
warmer future in high latitude sites [21] due to increasing soil frost 
depth associated with reduced snow cover or faster melt rates 
[22,23]. Winter soil frost could result in more fine root mortality 
[24] and intensify leakage of soil nutrients [25] . The increasing 
intensity of freeze-thaw cycles [26] could affect overall soil 
aggregate stabihty [27] and/or hydrological processes and flow 
paths as ice fiUs up the soil pore spaces [28] . The biogeochemical 
effects of intensified freeze-thaw cycles could include altered rates 
of microbial respiration, cell wall lyses and changes in the rate of 
decomposition of soil organic matter [29] . The latter could lead to 
increasing CO2 emissions from soils [13,30] resulting in a positive 
feedback to climate change [15]. The relative stability of carbon 
stored in high latitude forest catchments could be altered [31], thus 
shifting the boreal forest from sink to a source of carbon as 
temperature changes [32]. 

Therefore, there is an urgent need to properly represent soil 
temperature dynamics in forest carbon cycle models so as to 
constrain the uncertainty in future projections of environmental 
response to climate change. Here we present an extended version 
of the soil temperature model [33] used in the INC A suite of 
integrated catchment models [34—37]. The objective of this study 
is to address the question of how future climate change could affect 
soil temperature response in a set of well-monitored forest sites 
aligned in a south-north gradient in Sweden. Given the availability 
of long term soil temperature measurements at a series of depths in 
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the four Swedish Integrated Monitoring (IM) sites, we hope tliat 
this study will lead to improved soil temperature simulations when 
thermal conductivity is allowed to vary throughout the soil profile. 

Ethics Statement 

Soil temperature measurements were collected at the four 
Swedish Integrated Monitoring sites as part of the routine long- 
term monitoring conducted by the Swedish University of 
Agricultural Sciences, IVL, and others [38]. All sites have public 
access so no additional permissions were required. Our field 
studies did not involve endangered or protected species. All soil 
temperature data are publicly available through the Integrated 
Monitoring Program web site hosted by the Swedish University of 
Agricultural Sciences (http://webstar.vatten.slu.se/db.htnil) 

Study Area 

This study was conducted at the Swedish Integrated Monitoring 
sites (IM); Aneboda, Gardsjon, Kindia and Gammtratten 
(Figure 1). The Swedish IM sites are a set of well monitored 
headwater forest catchments and are part of a Europe-wide 
initiative to collect long time series of data on a number of key 
ecosystem variables. Data from the IM sites are used to study long 
term changes in fundamental catchment processes across the 
climatic gradient of Sweden. Due to the lack of land management 
such as agriculture or forestry, well-defined catchment boundaries 
and the availability of long time series of climate data, the IM 
catchments are suitable for the modeling exercise presented in this 
study. 

2.1 Specific catchment descriptions 

The Aneboda catchment (0.19 km^) is situated in the Sma- 
landska highlands (57°07'N, 14°03'E) and the biome is charac- 
terized as Boreo-nemoral (Figure 1). The mean air temperature in 
the catchment was 6.5°C and precipitation was 796 mm/yr 
(1996-2008). Norway spruce {Picea abies) is the dominant tree 
species (73%) with some Birch (Betula spp.) (20%). Other tree 
species with minor representation in the catchment include Scots 
Pine {Pinus sylvestris) (3%), Beech [Fagus sybatica) (1%) and Alder 
{Alms spp.) (2%). The soil is dominated by glacial tiU, with the 
proportion of wet soils amounting to 1 7 % on granite bedrock. In 
2005, the catchment was hit by a severe storm followed by a bark 
beetle infestation that caused extensive damage to the forest [38]. 

The Gardsjon catchment (58°03'N, 12°0rE) is located in 
Bohuslan in the Swedish west cost approximately 10 km from the 
sea (Figure 1). The catchment (0.04 km ) is the smallest of all the 
four IM catchments. The catchment biome is also characterized as 
Boreo-nemoral with mean annual air temperature of 7°C and 
annual precipitation of 1 1 1 1 mm/yr (1996-2008). The vegetation 
cover is dominated by Norway spruce (65%) but Birch (14%) and 
Scots Pine (17%) are also present. The catchment geology also 
consists of granitic bedrock underlying glacial till soils. Soils in 
Gardsjon are very shallow and interrupted by bedrock outcrops in 
some part of the catchment. The proportion of wet soils within the 
catchment is 10%. 

Kindia is situated towards the middle of Sweden (59°45'N, 
14°54'E) and has an area of about 0.2 km^ (Figure 1). The 
catchment biome is characterized as Southern-boreal with mean 
air temperature of 5.2°C and mean annual precipitation of 
854 mm/yr. The catchment also has significant bedrock outcrops 
(41%) and is quite steep with elevation ranging from 100- 
400 meters above sea level (m a.s.l). The vegetation is dominated 
by Norway spruce (83%) but Birch (14%) and Scots Pine (2%) are 
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Figure 1. Map of study catchments as aligned along the 
climatic gradient of Sweden. 

doi:10.1371/journal.pone.0093957.g001 

found in the catchment. The soil consists of 24% wet soil with mire 
in the center of the catchment. 

Gammtratten is the northernmost (63°5rN, 18°08'E) and 
largest (0.4 km^) of all the IM catchments (Figure 1) with mean air 
temperature of 2.4°C and annual precipitation of 579 mm/yr. 
The catchment is categorized as middle-boreal. Vegetation mainly 
consists of Norway spruce (70%), with some Birch (16%) and Scots 
Pine (13%). Pine is mosdy found in the higher elevations of the 
catchment, where there are several small mires. The percentage of 
wet soils was 16% with presence of granite bedrock underlying the 
glacial till soU. 

Climate 

3.1 Historical climate 

The data requirements for the soU temperature modelling and 
analysis presented here include observed average daily air and soU 
temperature. Continuous time series of both air temperature and 
precipitation (1996-2008) for each IM site were obtained from the 
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IM database at Swedish University of Agricultural Science. Future 
projections in Gammtratten were driven by climate data from the 
nearby Krycklan catchment [39]. Our assumption of using 
Krycklan was justified as both catchments have similar weather 
patterns (with of 0.73 for precipitation and 0.98 for air 
temperature). Due to large gaps in the soil temperature data, 
depths with the most consistent long term series were used for 
model calibration and validation in each catchment. Soil 
temperatures at three depths (Table 1) were modelled at each site. 

3.2 Future climate 

3.2.1 Climate models. Future climate data were based on 
ENSEMBLES project outputs [40]. The ENSEMBLES project 
utilized an ensemble of Regional Climate Models (RCMs) driven 
by different Global Climate Models (GCMs), to generate a matrix 
of possible climate projections under International Panel on 
Climate Change (IPCC) emission scenarios [1]. The ENSEM- 
BLES project generated a range of possible climate outcomes 
which help to make future projections more statistically reliable. 
Since the climatic factors (e.g. air temperature, winds, etc.) are a 
connected system in the atmosphere and extend all over the world, 
the GCMs operate at a global scale. The emission scenarios in the 
GCMs were based on assumptions about future population, 
economic development, etc. These assumptions were then 
translated into anthropogenic emission scenarios used as forcing 
for the GCMs, though different feedbacks forcing are also 
integrated include melting of ice caps [23] and changes in soil 
CO2 balances [12] amongst others. As a result, all GCMs are 
designed to model the complex climatological system of the earth 
with each GCM having certain processes and feedbacks repre- 
sented differently. This leads to differences between GCM 
representations of the future as a result of varied responses to 
forcing and differences in RCM responses even at much finer 
resolutions. 

3.2.2 Bias correction of climate models. Downscaling is 
often used in climate impact studies to create outputs with higher 
resolution [41—43]. Despite their finer resolution, RCMs are still 
coarse for direct applications in local impact studies as there are 
often biases between RCMs outputs and measured site specific 
conditions. Therefore we utilized distribution mapping in this 



Table 1. Catchments and their soil profile depths used In the 
modelling analysis presented in this study. 





Catchment 


Depth in soli profile (cm) 


Aneboda, top layer 


10 


Aneboda, middle layer 


32 


Aneboda, bottom layer 


58 


Gardsjon, top layer 


0* 


Gardsjon, middle layer 


10 


Gardsjon, bottom layer 


25 


KIndIa, top layer 


5 


Kindia, middle layer 


20 


Kindia, bottom layer 


35 


Gammtratten, top layer 


5 


Gammtratten, middle layer 


29 


Gammtratten, bottom layer 


40 



Study to correct biases in RCM-simulated temperature and 
precipitation. This approach has been shown to be the best bias 
correction method for small and meso-scale catchments in Sweden 
[44,45]. Here, 15 RCM ensembles were bias corrected to local IM 
conditions for the period 2061-2090 (Table 2). The general 
principle of distribution mapping is to fit cumulative distribution 
functions (CDFs) of observed data to CDFs of RCM outputs in the 
control period (1996-2008) according to Eqn. 1. 



'^coiitr ^ 



(F{Ccontr\pl contr contr 



where C is the climate variable of interest (precipitation or 
temperature), C* represents the bias-corrected climate variable, F 
stands for the theoretical CDF (Gamma or Gaussian) and for 
its inverse, pi and p2 are the distribution parameters (a and fi for 
Gamma distribution, fi and ct for Gaussian distribution), the 
subscripted expression ohs indicates observations and contr stands 
for the RCM-simulated control run (1996-2008). 

The derived distribution parameters from Eqn. 1 were then 
apphed to future series to adjust the RCM climate variables (Cjcen) 
according to Eqn 2. The subscripted expression seen indicates the 
RCM-simulated scenario run (2061-2090). The bias corrected 
ensemble RCM data were used as driving variables for the 
prediction of future soil temperatures at all depths in our study 
catchments except Gammtratten which was driven by bias 
corrected series from nearby Krycklan catchment. We refer 
readers to Teutschhein and Seihert [44] for more detailed descriptions 
of this technique. 



C* =F' 

seen 



' {F{Cscen \plconlr,p2conlr)\ploh.i,p2ohs), (2) 



Modelling Analyses 

4.1 Background on soil physics 

The soil temperature model was derived from two differential 
equations describing combined water and heat flows in seasonally 
frozen soil [46]. These equations were derived from the law of 
conservation of energy and mass (Eqns. 3 and 4) and the fact that 
water and heat spread out in the soil profiles along gradients of 
water potential and temperature (Darcy and Fourier laws). 



^ dT 81 8 

^'^-P''^'Jt=8z 



r(h\ ^'^ J. Pi ^ 



KriO) 



8T 

5^ 



8h 



- Cwq\ 



8T 



-S(h) 



(3) 



(4) 



* Simulated as 1 cm depth. 
doi:10.1371/journal.pone.0093957.t001 



where Z (m) is a space coordinate, T (°C) is soil temperature, Kx 
(Wm ■ °C ) is the soil thermal conductivity, Z/r (f kg ) is latent 
heat of fusion and water, Q [f m ■'*°C ') is volumetric specific 
heat of the soil, p/ (kg m '^) is the density of ice, pjy (kg m ^) is the 
density of water, g^/ (m s ') is flow of water, q (dimensionless) is a 
volumetric water content, / (dimensionless) is a volumetric ice 
content, h (m) is the soil water potential, C(h) (m ) is the 
differential moisture capacity, K(h) (m s ') is the unsaturated 
hydraulic conductivity of the soil matrix and S(h) (dimensionless) 
represents the water taken up by the roots. These equations 
calculate water and heat flow based on soil properties; the water 
retention curve, unsaturated and saturated hydrauhc conductivity. 



PLOS ONE I www.plosone.org 



3 



April 2014 I Volume 9 | Issue 4 | e93957 



Soil Temperature Dynamics across Sweden 



Table 2. Underlying RCMs for the bias corrected site-specific climate scenarios, their notation, resolutions, driving GCMs, emission 
scenarios and performing institutes. 



No. 


Notation 


Institute 


RCM 


Resolution 


Driving GCM 


Emission scenario 


1 


C41 HAD 


C4I 


RCA3 


25 km 


HaHCM^nifi 
nduv^ivi jv.^ 1 o 


Al B 


2 


DMI ARP 


DMI 




25 km 


ARPEGE 


Al B 


3 


L.'IVII DAIVI 


DMI 




25 km 


BCM 


Al B 


4 


DMI Fru 

L/FVLI tv_n 


DMI 


HIRHAMS 


25 km 


FCHAMS 


Al B 




ETHZ 


ETHZ 


CLM 


25 km 




Al B 




ur HAnn 


HC 


ridLinivi jvju 


25 km 




Al B 


7 


ur UIAD^ 
rn^ riML/ J 


HC 


ridLinivi 1 D 


25 km 


ndUv.ivi3v^ 1 D 


Al B 


8 


HC_HAD16 


HC 


HaclRM3Q3 


25 km 


HadCM3Q3 


AlB 


9 


KNMI 


KNMI 


RACMO 


25 km 


ECHAM5 


AlB 


10 


MPI 


MPI 


REMO 


25 km 


ECHAM5 


AlB 


11 


SMHLBCM 


SMHI 


RCA 


25 km 


BCM 


AlB 


12 


SMHLECH 


SMHI 


RCA 


25 km 


ECHAM5 


AlB 


13 


SMHLHAD 


SMHI 


RCA 


25 km 


HadCM3Q3 


AlB 


14 


CNRM 


CNRM 


Aladin 


25 km 


ARPEGE 


AlB 


15 


ICTP 


ICTP 


RegCM 


25 km 


ECHAM5 


AlB 



Note that both the numbering sequence and/notations are used interchangeably to refer to each RCM throughout this study. 
doi:l 0.1 371 /journal.pone.0093957.t002 



heat capacity (including latent part during thawing/melting) and 
thermal conductivity [33]. 

4.2 Rankinen model and the extended version. The 

model used in this study was based on the soil temperature model 
described in Rankinen et al. [33] which was based on the 
simplifications of Eqns. (3) and (4). The model calculates soil 
temperature at different depths at a daily time step with the fuU 
consideration of the influence of snow cover. Simplifications made 
by Rankinen et al. [33] neglected the influence of changes in soil 
water content on soil temperature. This made it easier to solve 
Eqn. (4) and aU water related terms from Eqn. (3) with the 
exception of the ice term. 

In evaluating the derivatives from Eqn. (3), boundary conditions 
were set as Tguj^p (temperature at the surface, replaced by air 
temperature T^jr in Eqn. (5), T7 (temperature in the soil evaluated 
according to space reference ^ and T; which is soil temperature at 
2^. For purposes of simplification, this last term was assumed to 
equal T^, indicating that heat flow under the soil layer of interest 
could be taken into consideration. The derivative of the soil 
thermal conductivity in the soil profile was replaced by a constant, 
representing the average thermal conductivity of the soil. 
However, the soil thermal conductivity was linked with the soil 
moisture [46], which varies throughout the soil profile and 
throughout the seasons. As a result of this simplification, the model 
lost its validity under very wet or dry conditions but gready 
simplified estimation of soil temperature (Eqn. 5). 



Tl 



(Cs + C/c£)(2Z, 



^2 ["^'aIR ^z] 



(5) 



As eqn. (5) did not take the influence of snow into account, the 
equation was extended by an empirical relationship (Eqn. 6) that 
could simulate the insulating the effect of snow cover and soil 
temperature could thus be calculated (Eqn. 7). 



r'+l . 



-T' 



-fsDs 



(6) 



r'+l - 



■T'+- 



{Cs + C,ce)(2Zs 



o-fsDs^ 



(7) 



where Kr 

-3. 



(W'm"'°C"') is the thermal conductivity, Cs 
(J m~''°C~') is the specific heat capacity of the soil, Ciqe 
(J m """C ') is the specific heat capacity due to freezing and 
thawing (latent part) as well as _/J (m~') which represents an 
empirical snow parameter. These represent model parameters that 
can be used in the model calibration. 7^/^ (°C) and Ds (mm) 
represent air temperature and snow depth. 

This study further explored the possibilities of improving the soil 
temperature simulations (Eqn. 8) when the thermal conductivity is 
allowed to vary throughout the profJe. This extended version also 
utilized temperature inputs below the soil layer of consideration. 
Therefore, parameters controlling the lower soil thermal conduc- 
tivity Ky_low> lower soil specific heat capacity C\low, and lower 
soil temperature Tixiw'^^^^ also added to the model (Eqn. 8). 



r'+l . 



AtKrj 
Cs,i2{Z, + zf 



■T'z] 



(8) 



T,-T' 



where is a constant indicating the space coordinate for the lower 
temperature influence. C = Cs, T'^>0 or C=Cs + Cice', 7z < 0 
The model was implemented in MATLAB. 

4.3 Model simulations and projections 

The model was calibrated to present day condition (1996- 
2008). Since there were gaps in soil temperature observations, the 
calibration and validation periods were not split by a specific date 
but by examining the number of observations available and then 
splitting the data series in thirds. The first two thirds of the series 
were used for calibration, making the calibration and validation 
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Aneboda 10 cm 




Gammtratten 5 cm 




Aneboda 32 cm 



Gammtratten 29 cm 





Aneboda 58 cm 



Gammtratten 40 cm 





" Observed 



" Modelled 



" Treshhold for calibration 



" Observed 



" Modelled 



" Treshhold for calibration 



Figure 2. Simulated and observed soil temperatures for the representative southern (Aneboda) and northernmost (Gammtratten) 
catchments. 

doi:1 0.1 371 /journal.pone.0093957.g002 



periods difTer between catchments and soil depths (e.g. Figure 2). A 
Monte Carlo sampling technique was used for model optimization 
by random sampling of parameter spaces as well as tracing out the 
structure of the model output [47] . This allows proper evaluation 
of multiple combinations of parameter settings in models with 
significant uncertainty in inputs and outputs. By changing one 
parameter at a time (as is usually done in manual calibration), 
model performance might be less than optimal as there would only 
be one degree of freedom. However, evaluating multiple 
combinations of parameters manually is often impossible in more 
complex systems. 

Parameters that have a physical interpretation (grey box 
modeling) must be combined with prior knowledge of the system. 
In the soU temperature model described in this study, six out of the 
seven parameters (Table 3) have clear physical interpretations as 



their values are limited by their physical boundaries. It is therefore 
of importance that the model parameters are set within a 
reasonable range when performing Monte Carlo analysis 
(Table 3). Choosing the parameter range for individual sites and 
depths was difficult. However, there are a number of strategies in 
choosing a suitable parameter ranges. These include literature 
values, experimental results and measurements as well as expert 
judgments. The range at which the upper soil temperature 
parameters (C?, K'j-,fs and Cjce) were allowed to vary was based on 
the parameter range proposed by [33,48]. The range of the lower 
soil temperature parameters {Ti^ow, Cs,LOW and Ky^low) were 
vmknown, and were therefore set according to judgments. 

The Monte Carlo analysis was conducted by running the model 
100 000 times, sampling a new set of randomized parameters for 
each model run from the chosen parameter ranges. This process 



Table 3. Parameter ranges for the Monte Carlo simulations. 



Parameter 


Unit 


Monte Carlo ranges 


Specific heat capacity of soil, Cj 


J m C 


0.5-3.5 (10"^) 


Soil thermal conductivity, 


Wm"'°C"' 


0-1 


Specific heat capacity due to freezing and thawing, C/ce 


1 ™-3c, 1 

J m C 


4-15 (10*^) 


Empirical snow parameter, fs 


m 


0-10 


Lower temperature, T^qw 


°C 


0-1 


Thermal conductivity, lower part, Kjiow 


Wm"'°C"' 


0-1 


Specific heat capacity of soil, lower part, Cs,low 


J m C 


0.5-3.5 (10"^) 



doi:1 0.1 371 /journal.pone.0093957.t003 
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Aneboda 
Gardsjon 
Kindia 

Gamnntratten 



X 10 



1.6 

7 



Figure 3. Cumulative relative distributions for Kr{p\ot a) and Q/cf (plot b) from the Monte Carlo simulations of middle soil layers 
across the four IM catchments. 

doi:1 0.1 371 /journal.pone.0093957.g003 



was repeated for each site and for all depths. The high number of 
simulations increased the likelihood that the whole parameter 
range was well sampled. Each parameter was sampled using the 
MATLAB randQ command, scaled to fit each parameter range. 
Model performances during calibration period were measured 



using the Nash and SutclifFe (NS) statistic [49] while validation 
period performance was measured using a combination of NS, 
(calculated based on linear fitting of observed and modeled values) 
and root mean square error (RSME). The Monte Carlo outputs 
were further analyzed using cumulative distribution frequency 




Figure 4. Ensemble projection of air versus soil temperature across the soil profile in the four integrated Monitoring catchments. 

doi:10.1371/journal.pone.0093957.g004 
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Figure 5. Corresponding change in tKie future air and soil temperature relative to the current conditions (shown in Figure 4). 

doi:1 0.1 371 /journal.pone.0093957.g005 



(CDF) of the top 5000 simulations. A sensitive parameter tends to 
have a non-uniform posterior distribution while a non-sensitive 
parameter distribution spreads over the entire range. Parameter 
values for the behavioral runs (top 100) and bias-corrected climate 
data were used for projecting future conditions (Figure 3). The 
projected series were further analyzed on both annual and 
seasonal scales for each of the 15 RCMs during 2061-2090. 

Results 

5.1 Soil temperature simulations 

The model successfully captured inter-annual variation in soil 
temperature at all catchments. Figure 2 shows the model 
calibration in Aneboda and Gammtratten to illustrate the south- 
north gradient. Model performance in Aneboda ranges from NS 
value of 0.96-0.97 (calibration) and 0.96-0.97 (validation). The 
model also performed well in Gardsjon (NS 0.97-0.98), Kindla 
(NS 0.95-0.97) and Gammtratten (NS values of 0.96-0.98). There 
was greater variability in upper and middle layer soil temperatures 
at Aneboda than at Gardsjon. The variability increased northward 
and was consistent in the upper soil layer with the increasing 
influence of air temperature. The model failed to capture some of 
the winter soil temperature dynamics (Figure 2) at Gammtratten, 
making this catchment a unique site with clear insulating effects of 
snow cover as the soil temperature flattens out during winter 
months (Figure 2). 



5.2 Uncertainty analysis 

The vmcertainty analysis revealed parameter equifmality. This 
means that there is no single best model parameter set and that 
many model state descriptions can generate equally good 
calibration outputs [50]. For example, parameter CsjcE was 
highly insensitive (Figure 3), making it hard to identify optimum 
values of this parameter to best represent present day soil 
conditions. Additionally, it was hard to make cross-catchment 
comparisons since soil temperatures were simulated at different 
depths for the diflFerent sites. However, K7- was clearly the most 
sensitive parameter that differentiates each catchment along a 
south-north gradient in Sweden (Figure 3). Both Kindla and 
Gardsjon have narrower K-p ranges compared to Aneboda and 
Gammtratten. This shows the importance of soil thermal 
conductivity in regulating soil temperature in high latitude 
catchments. 

5.3 Ensemble projections 

5.3.1 Changes in air temperature. The RCM ensemble 
projected a range of possible future climates across Sweden 
(Figure 4). Changes in air temperature are projected to be greatest 
in the colder months (November-April) with the magnitude of 
change increasing northward. Future air temperatures <0°C were 
projected to decrease markedly in the northern catchments with 
possibility of complete disappearance of winter air temperatures 
<0°C in the southern catchments (Figure 4). The largest change in 
air temperature would be observed in winter (about 3°C rise) with 
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Figure 6. Mean annual mean air-soil temperature for the control period (calibration/validation period) using median year (red 
diamond) versus mean annual air-soil temperature ensembles (plain diamonds. Label 1-15 denote each RCM and fully described in Table 2), 
using ensemble median for Aneboda (A), Gardsjbn (B), Kindia (C) and Gammtratten (D) middle soil layer. Error bars represent corresponding annual 
standard deviations. 
doi:1 0.1 371 /journal.pone.0093957.g006 



a possibility of autumn warming up to 2.5°C (Figure 5). However, 
Gardsjon showed the possibility of sporadic warming up to 3.4°C 
in August, thereby shifting the warmest month of the year away 
from July in that catchment. The patterns of future change in air 
temperature in Kindia are fairly similar to the patterns projected 
in Gardsjon (Figure 5). There could also be an abrupt rise in air 
temperature up to 3°C in the month of August with the biggest 
winter air temperature change 3°C. 

The northernmost catchment showed markedly different 
magnitude and patterns of change compared to the rest of the 
southern catchments. Ensemble projections showed that winter air 
temperature change could be largest (3.4°C rise) in Gammtratten 
region with a possibility of sporadic warming events up to 3.9°C in 
February and May. Autumn warming, though slightly visible in 
other catchments, could become more common in this region in 
the future (up to 3°C rise). However, there could also be sporadic 
warming events peaking in November (3.9°C rise) before the onset 
of the colder winter months (Figure 5). These results show the 
possibility of increased variability and fluctuations of future climate 
across Sweden. 

5.3.2 Changes in soil temperature. The overall ensemble 
showed markedly different winter soil temperature responses to 
future climate change across the soil profile in each catchment 



(Figure 4). Though soil temperature responses increase northward, 
the pattern of change varies sigiiificantiy across months in each 
catchment (Figure 5). Despite their southern location, patterns of 
soil temperature response in Gardsjon are different from Aneboda. 
Changes in soil temperatures in Gardsjon were less variable and 
more synchronous with air temperature (Figure 5). Although air 
temperature could also be changing more rapidly than soil 
temperatures in this catchment on an annual scale (Figure 6), the 
largest change was projected to occur in the top soil with 
differences between the top/lower layers responses being most 
pronounced in late summer and early autumn. 

While the two southern catchments had relatively moderate soil 
temperatures throughout the year, changes could be more 
pronounced in Kindia and Gammtratten sites, the more northern 
sites (Figure 5). The biggest soil temperature changes in Kindia 
were projected to occur in April-August and the effect was most 
pronounced in the topmost layer relative to the middle and bottom 
layers which showed fairly similar response patterns. Middle and 
bottom soil temperatures could rise by about 1.5°C (about 2.5°C 
for topmost layer) in the future for most months of the year, except 
December-February where the model projected a temperature of 
about 0.5°C. 
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Despite the largest air temperature changes in the northernmost 
Gammtratten catchment, there was no corresponding soil 
temperature responses in the winter unlike other catchments 
where soil could get warmer during the cold and dormant winter 
season. This clearly shows the insulating effects of extensive snow 
cover in this region. In contrast to the almost unchanging in soil 
temperature in winter months, there was a steady rise in soil 
temperature after March and this could reach up to 4°C in May in 
Gammtratten (the largest change in all catchments). The 
subsequent stabilization in the summer months (about 2.5°C rise) 
is similar to the magnitude of changes observed in the other 
catchments. This suggests that soil warming could be largest in the 
spring season in die northern boreal region. 



climate impacts, they might be uninformative in high latitude 
boreal catchments unless complemented by seasonal projections. 
Our seasonal analysis showed that impacts could be more 
pronounced during winter months and the magnitude of the 
change in air temperature increases northward. This is because the 
period with winter air temperatures <0°C was projected to 
shorten considerably in the northern catchments that are 
characterized by significant snow cover which may disappear 
completely in the south. This will have large impacts on 
ecohydrology and biogeochemical processes including snow 
accumulation and melting, soil frost [53], ecosystem and forest 
productivity [19], changes in water quality [5], distribution of tree 
species and community composition [54] amongst others. 



Discussion 

6.1 Climate variability and change 

Several long-term monitoring studies have shown a possible 
warmer world toward the end of 20th century [3,5,9] . This trend 
might continue in the future if the current levels of human 
activities persist. The modelling analyses presented in this study 
showed that ensemble projections performed better in constraining 
the uncertainty in future climate change impacts across the 
climatic gradient of Sweden. This is because no single model can 
accurately represent the magnitude of change expected in the 
fiiture considering the level of uncertainty inherent in our current 
forecasting capability. However, the sign of the changes observed 
from this study and by others [20,51] could give some insights on 
how to devise plausible mitigation strategies for an uncertain 
future. Although our projected changes on annual scale are large, 
they fall within the range of values projected by other studies. For 
example, Kellomdki et al. [52] projected an up to 4.5°C rise in air 
temperature over Finland amidst other projections in Sweden 
[44,51]. While annual averages provide a quick overview of 



6.2 Soil temperature modelling 

The extended soil temperature model presented in this study 
offers improved soil temperature simulations as heat-flow are 
allowed to influence the soil layer of consideration from below. 
This is very important in understanding other biogeochemical 
processes as soil temperatures are becoming represented in other 
environmental source assessment models [36,42,55]. Recently, 
Winterdahl et al. [56] observed improved simulations of stream 
dissolved organic carbon concentration when soil temperature was 
integrated to the riparian processes. This showed a need for proper 
accounting for soil temperature in biogeochemical studies [5 7] and 
to further constrain uncertainty in biogeochemical models [34,35] 
for credible assessments of possible futures. 

The deepest soil layers consistendy showed less variability than 
the topmost layers. This is explainable as topmost layer can be 
influenced by fluctuations in daily air temperature. This is 
particularly notable in Kindla and Gardsjon where bedrock 
outcrops [38] can conduct additional heat from the atmosphere to 
the soil. In general, Gardsjon has shallower soils that can make the 
top two layers more sensitive to daily air temperature variability. 
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Figure 7. Summary of projected change in long term annual soil temperature using the RCM ensemble median (white diamonds) for 
Aneboda (A), Gardsjon (d), Kindla (K) and Gammtratten (C(\ versus corresponding observed median soil temperatures for the test 
period (red diamonds). Error bars represent corresponding annual standard deviations. 
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However, the large variation in the earlier observations at 
Aneboda could not be explained in terms of soil physical 
properties related to the proportion of soil wetness (17% in 
Aneboda; 24% in Kindla). Southern Sweden was hit by a severe 
storm in 2005 followed by pest infestation [38] that destroyed 
some forest vegetation (more severe in Aneboda than Gardsjon). 
However, no generahzation could be made on the impact of the 
storm in Aneboda or influence of soil heterogeneity in each 
catchment as soU temperature measurements were not made 
throughout the landscapes. 

The soil temperature model presented here successfully 
captured the inter-annual and seasonal dynamics of soil temper- 
ature in each catchment. The high NS values for model calibration 
are explainable as soil temperature is less variable unlike air 
temperature noted above. Though high summer and low winter 
soil temperatures were well simulated, the model slightly 
underestimated soil temperature in some years especially in 
northernmost catchment. This was particularly notable during 
autumn cooling and spring warming. Gammtratten was also the 
study site where the influence of snow cover change was clearly 
visible. This is not surprising in this region as similar observations 
were also reported in neighboring Finnish catchments [33]. 

General uncertainty anah sis showed that the model showed 
equifinality [50] due to strong correlation between parameters. 
However, the parameter Csjce appeared to be highly insensitive. 
This observation is not surprising since CsjcE only comes into play 
when soil temperatures drop below zero. This argument also holds 
for parameter fs, which only affects NS when there is snow. 
Calibrating only for the winter would reveal more about the true 
sensitivity of the mentioned parameters and this approach would 
be attempted in future work. However, sensitivity of Kx showed 
the importance of soil thermal conductivity and that an 
assumption of uniformity across soil profiles or catchments should 
not be made. This is important as Euskirchen et al. [11] clearly 
showed the significance of change in soil thermal conductivity on 
growing season length and forest productivity across boreal 
ecosystems. 

6.3 Ensemble projections 

This study shows that the relationship between air and soil 
temperature responses is not linear and that complex, climate- 
related feedbacks will increase in northern boreal biomes. Our 
projected range of future soil temperature responses is consistent 
with other studies including a recent study [58] which projected an 
up to 3.3°C rise in soil temperature in southern Quebec forest 
catchments in Canada, ^hang et al. [9] previously used the process- 
based NEST model across a climatic gradient in Canada and 
showed that soil temperature (at 20 cm depth) increased by 0.6°C 
for 1°C rise in air temperature. They attributed the difference to 
the effect of more insulation from snow cover in winter in high 
latitude catchments. Our study also indicated that air temperature 
could increase more rapidly than soil temperature across the 
climatic gradient of Sweden. This could make future soil/air 
temperature in northernmost Gammtratten shift towards present 
day conditions in the south by the end of the century (Figure 7). 
One possible implication of the results presented here is a shift in 
biomes from middle-boreal that characterize the northernmost 
Gammtratten catchment to southern- boreal or in the worst case, 
boreo-nemoral that characterize southern Sweden. 

The large soil temperature changes and possibility of lesser 
snowfaU in the southern catchments due to warmer winter will 
have biogeochemical implications for peat)^ soils. Increased winter 
soil wetness and less frost could reduce the carrying capacity of the 
soil for winter forestry operations [52]. Though drier soils could 



support deep winter soil frost under warmer conditions [26,53] but 
possible persistent soil frost in the north cottld also be of ecological 
concern. For example, our results showed possibility of minimal 
winter soil temperature changes in Gammtratten catchment. 
Although this remains around the freezing point (Figure 5), there is 
a possibility of slightly colder soils. This borderline winter soil 
temperature responses around the freezing point in Gammtratten 
is not surprising as Campbell et al. [59] recently reported similar 
results from the Hubbard Brook Experimental Forest in the 
United States. Their SHAW model similarly showed some biases 
in partitioning precipitation to either rainfall or snowfall around 
the freezing point. They suggested that precipitation would be a 
less important factor influencing snowpack and soil frost in the 
winter in the future. This partly justifies our lack of snowfall data 
to constrain this uncertainty in our model but the tendency of 
decreasing winter soil temperature at Gammtratten is not 
unambiguous as whether it would result from less or more snow 
cover. Several studies have also shown possibility of colder winter 
soils as the future gets warmer [22,3]. Therefore, warmer air 
temperature can be a limiting factor for deep soil frost as length of 
frost days cotdd be concurrently reduced. Decreasing frost days 
are evident at Gammtratten where the model showed possible 
largest increase in soil temperature in the spring. This is an 
indication of a faster snowmelt regime and as a result could lead to 
extended growing season in the future. Recent analysis of 
historical data by [5] had shown already increasing trends in the 
growing season length in this region 

Stromgren and Linder [19] conducted a soil warming experiment 
and showed that heated plots also had less snow cover and favored 
more stem volume production. This was due to possible extension 
of growing season [1 1] and ambient environment that favors more 
nutrient availability [14]. Bergh and Linder [60] also conducted a 
proxy study of future climate effects by burying heating cables in 
the soil. They showed that snowmelt and soil thawing occurred 
earlier in heated than unheated plots. Several other snow 
manipulation experiments have also shown that the persistent 
deep soil frost could increase organic matter leaching [10] and 
lower rates of winter respiration [13]. A drastic reduction in 
ceUtjlose decomposition (up to 46%) was observed in snow 
removal experiment in northern boreal forest in Sweden [29]. 

To further test the internal working proct'ss of the model 
presented in this study, application of the model to catchments 
with both soil temperature and snow data would be worthwhile to 
constrain the uncertainty in model structure and future. Compar- 
ing the soil temperature measurements from the same depth in the 
soil profiles could give more credible comparison between sites to 
evaluate whether parameter behaviors are due to soil type, or 
changing soil properties down the profile. While this is important 
to constrain the model uncertainty, it might be difficult to attain. 
One recent study [61] showed that soil properties can significantiy 
affect thermal distribution in soil profiles such that direct 
comparison could be difficult even if the measurements were 
taken at the same depth in each catchment. Future work should 
also consider using only winter soil/ air temperatures data for 
parameter estimation when simulating winter temperatures. The 
overall conclusion of this study showed that the effects of climate 
change on soU temperature responses in snow dominated region is 
complex. General assumptions of soil temperature dynamics based 
on future air temperature change alone are not enough and should 
be avoided in high latitude ecosystems. 
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